A Particle-Particle, Particle-Density (P^D) algorithm for the 
calculation of electrostatic interactions of particles with slab-like 

geometry 

S. Alireza Ghasemij*] Alexey Neelov, and Stefan GoedeckeJl] 

Condensed Matter Theory Group, Department of Physics & Astronomy, 
Klingelbergstrasse 82, Basel 4056, Switzerland 
(Dated: February 1, 2008) 

Abstract 

We present a fast and accurate method to calculate the electrostatic energy and forces of interact- 
ing particles with the boundary conditions appropriate to surfaces, i.e periodic in the two directions 
parallel to the surface and free in the perpendicular direction. In the spirit of the Ewald method 
the problem is divided into a short range and long range part. The charge density responsible for 
the long range part is represented by plane waves in the periodic directions and by finite elements 
in the non-periodic direction. Our method has computational complexity of 0{Ng \o^{Ng)) with a 
very small prefactor, where Ng is the number of grid points. 
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I. INTRODUCTION: 



Simulations of systems with slab-like geometries are of great importance. Problems in- 
volving surfaces, interfaces, tip-surface interaction in scanning probe microscopy simulations, 
electrolytes trapped between two plates, thin films of ferrofluids, etc. all fall into this cate- 
gory. Calculating the Coulomb interactions in such setting is a major challenge. With free 
boundary condition (i.e. the potential tends to zero at infinity) the scaling of the trivial 
direct summation is 0{N'^) where is the number of particles. In the case of 2D periodic 
and ID free (2DP1DF) boundary conditions (BC) the situation is even worse. In principle 
one would then have to include into the summation the interations with all the periodic 
images in the two periodic directions. 

P n 

Algorithms such as Ewald-based methods pj], fast multipole methods(FMM)[2], P^M 
method .3], and convergence factor approachesjj, [s, O] have therefore been generalized to 
2DP1DF problems. Handling different types of BC in EMmI?] is straightforward. In addition 
the FMM methods have the ideal linear scaling. Unfortunately the prefactors in FFM 
methods are typically large and so the FMM methods are in many cases only faster than 
other methods for > 10^, where N is the number of particles. Another drawback of 
FMM that is important in molecular dynamics is that the approximate FMM forces are 
not analytical derivatives of the approximate energy. Therefore the energy is not conserved 
during the molecular dynamics simulation. High accuracy energy conservation is therefore 
impossible. 

Ewald methods for 2DP1DF boundary conditions, called EW2D, have been developed 



Refs. |8|, y, b^. A comparison of three versions of EW2D method can be found in 
Ref. .Unfortunately, the practical use of the EW2D sum is hampered by the occurrence 
of a reciprocal space term. The resulting Fourier space sum does not allow for a product de- 
composition as it is done in the three-dimensional periodic Ewald method and therefore the 
method has a scaling of 0{N'^). In 2002 Arnold and Holm developed MMM2Dll2|(MMM 
with 2DP1DF BC), which is found to be the best in terms of accuracy. Another advantage 
of this method is that it has "a priori^ error estimates. However, because of its 0{N3') 
scaling it is only suitable for a small number of atoms. 

A rather simple approach is to use the standard three dimensional periodic Ewald method 
(EW3D) also for 2DP1DF boundary conditions. Spohr showed that the regular EW3D 
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method almost reproduces the EW2D results [isl. provided that the box length in the non- 
periodic direction is about five times larger than those in the periodic directions and that 
there is empty space of sufficient thickness in the basic periodic box to dampen out the inter- 
slab interactions. There are also methods with correction terms to make the 3D periodical 
scheme applicable to the 2DP1DF systems and resolve the problem of slow convergence with 
respect to thickness, so that a medium size gap(empty space) is enough. The EW3Dc[l^,Q 
method consists of a modification of EW3D to account for the slab geometry and addition of a 
correction term to remove the forces due to the net dipole of the periodically repeating slabs. 
Methods with layer correction terms to eliminate the inter-slab interaction, in addition to the 



correction term responsible for net dipole, have been mixed wit 



16|, 



17j, P^MLC 



1 mesh-based methods, thus 
16,13. The main drawback 



almost linear scaling is achieved e.g. EW3DLC 
of these methods is that the errors in the forces on the particles near to the surfaces are 
more than in the middle. 

In this paper we present a method which fills the gap of absence of an efficient method 
for medium size systems having 10^ — 10^ particles. Because our method is not based 
on a modification of a fully periodic method, no replication is needed in the non-periodic 
direction, leading to smaller memory and CPU requirements. In contrast to some others, our 
method does not impose any restriction on the distribution of particles in the non-periodic 
direction. 



II. COULOMB INTERACTION FOR SYSTEMS WITH 2DP1DF BC 

Consider a system of N particles with charges qi at positions in an overall neutral and 
rectangular simulation box of dimensions L^iLy and L^. The Coulomb potential energy of 
this system with periodic boundary condition in two directions and free boundary conditions 
in the third direction(let us say in the z direction) can be written as 

/ N 

n i,j=l ■' 

where Vij = r j — rj and n = {n^Lx, nyLy, 0), with Ux, Uy being integers. The prime on the 
outer sum denotes that for n = the term i = j has to be omitted. 

In the Ewald-type methods the above very slowly converging sum over the Coulomb po- 
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tential function is split into two sums which converge exponentially fast, one in real space 
and the other in the Fourier space. This splitting can be done by adding and subtracting 
a term corresponding to the electrostatic energy of a system of smooth spherical charge 
densities,pi(r), centered on the particle positions: 



/ N 



E 



1 ^' 

4ee 



\Vij + n| 



p,(r)pj-(r' + n) 



dvdr' 



n i,j=l 
N 



Pi{r)pj{r' + n) 



r — r 



dvdv' 



r — r 



i=l 



(2) 



The aim of the last term is to subtract the self energy for n = and i = j which is included 
in the second term. 

Even though Ewald-type methods allow for any choice of Pi(r), it was noted in Refs. HQ 
that Gaussians are virtually optimal in practice. Choosing Pi(r) to be a Gaussian function 



— exp 



r - Yi 



a' 



(3) 



leads to a well-known formula for the first and the third term in Eq. ([2]) . 



E 



I N 

EE 

n i,j=l 
N 





erfc 


>ij+nr 






rij + n 





4ee 

n «j=l 
1 ^ 



+ 



Pi(r)pj(r' + n) 



r — r 



dvdv' 



av 27r 



i=l 



(4) 



Obviously, the calculation of the third term is trivial. Since the interaction in the first 
term is decaying exponentially it can be made of finite range by introducing a cut-off. The 
error resulting from the cut-off is then also exponentially small and the short range term 
can be calculated with linear scaling. We have calculated the short range part and also the 
contribution of forces from long range as it is described in Ref. 18 1 
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The major difficulty is the calculation of the second term. A method to solve the Pois- 
son's equation under 2DP1DF boundary conditions has recently been put forward by L. 



Genovese |20|. Our approach is similar. As in Ref |20| we use plane waves [2 11] to represent 
the charge density in the periodic directions. Whereas Genovese et al used scaling func- 
tions as the basis in the non-periodic direction, we use finite elements for that purpose. 
Scaling functions are presumably the optimal choice in the context of electronic structure 
calculations where the charge density is given on a numerical grid. In our case the charge 
distribution is a sum over smooth Gaussians that can easily be represented by our mixed 
basis set of plane waves and finite elements. As will be seen we can avoid storing any kernel 
if we solve a differential equation along the z-axis instead of solving an integral equation. 
We will use a family of finite elements that allows to solve the linear system of equations 
resulting from the differential equation very efficiently. 

A. Calculating the long range part 

The second term in Eq.Q, can be written as 



Eion, = ]: [ p(^)(r)\/(r)rfr (5) 



where 



N 



p(^)(r) := 5^p,(r) (6a) 



1=1 



V(r) := I /^dr' (6b) 
i5R3 |r - r'l 



N 



n j=l 

We consider a system with a charge density that is only localized in the non-periodic di- 
rection, in our notation z; p{x,y,z) = \/{x,y,z) E \ z ^ [zi,Zu]. We define the cell 
containing the continuous charge density as: 

V := [0,LJ ® [0,Ly] [zi,Zu] 
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In our case the length of V in 2; direction Zu — zi is plus twice the cut-off for Gaussians. 
Since p(r) is periodic in x and y direction, V{r) is periodic too, so we can rewrite Eq.([5]) as: 



long 



p{r)V{r)dr 



(7) 



and V{r) can be calculated in an alternative way to Eq. fl6b|) . It can be considered as 
the solution of Poisson's equation with 2DP1DF BC: 



VV(r) = -47rp(r) 



(8) 



The charge density and the potential are periodic in x and y directions. Hence we can 
write the potential and the charge density in terms of Fourier series: 



V{x,y,z)= ^ Ckiiz)exp 

k,l=—oo 

p[x,y,z)= -ZJ—^^P 



k,l=— 00 



2tn{— + 

,kx ly 
2m(— + / 



(9a) 

(9b) 



Inserting Eqs. (pa|) and (l9bl) in Eq.® yields: 



d2 



ill] Cki{z) = riki{z) 



(10) 



f^2 p 

Jki:=2nJ— + — 



Vkliz) 



-Air 
LxLy Jo 



X exp 



p{x,y,z) 



-2mi— + 



dxdy 



(11) 



To solve the differential equation (fTOj) one needs to have boundary conditions at z ^ ±00 
for Cki{z). The potential obtained by solving Poisson's equation should be the same as the 
one in Eq. (16bp . Hence we derive the boundary condition in the non-periodic direction from 



Eq. (I6bl) . By performing the Taylor expansion of pc-jy about 2;' = in the integral expression 



of Eq. (l6bl) for the exact potential V{x,y,z) arising from our periodic charge distribution 
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p(r) 



Zu poo POO 



V{x, y,z) = I I I dx'dy'dz' —; — 



k,l=—oc 



Zl J — oo J — oo 



E l kiy^ ; 



Lx Ly 



(12) 



one can show that V{x, y, z ^ ±00) = =p/3 where (3 is proportional to the dipole moment of 
the charge distribution along the z direction 

P=l r\oo{z')z'dz' (13) 



For the Gaussian charge distributions given by Eq. the above integral can be calculated 
analytically and f3 is calculated exactly. 



-2vr ^ 



P = -r-rT.^^'^ (14) 



x-^y 



i=l 



This boundary condition for the potential gives the following conditions for the 7's. 

• 7 = 7oo = ^ ■^coo{z) = rioo{z) We solve this differential equation with boundary 
condition Cqq^z — » ±00) = 

• 7 = 7a:/ 7^ ^ — 'jIi^ Cki{z) = rjki{z) For all of these differential equations we 
have to impose BC of the form Cki{z ±00) = 0. 

The solution for 000(2:) is a linear function outside the interval [zi^Zy\. Since the boundary 
conditions are applied at infinity the linear term has to vanish and one has to satisfy Dirichlet 
BC for Coo, namely Cqq{zu) = — /3 and 000(2;;) = p. For \k\ + \l\ > 0, Cki{z) will have Robin 
BC as explained below. The potential is thus not modified if one takes for instance a 
computational box that is thicker in the z direction than necessary. The thinnest possible 
box is the one that just includes the region where the charge is nonzero. 
For z G {—Qc, Zl] we have rj^i^z) = thus it yields 

c{z) = c{zi)e""^'^'"''^ (15) 
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Both c{z) and its derivative must be continuous. So performing left differentiation at zi we 
get: 

c'{zi) - 7kic{zi) = (16) 
With a similar procedure we obtain the BC at Zu'- 

c{zu) + 'jkic{zu) = (17) 
These BCs are in agreement with the BCs resulting from the Green functions in Ref. ^ 

B. Solving the ordinary differential equation using the finite element method 

We recapitulate the procedure of solving the differential equation for the case + |/| > 0, 
i.e. 7fci 7^ 0, using the finite element method. For the case k = I = the approach is similar, 
with the only difference that the Dirichlet BC are used. The case k = I = can be found in 



many manuscripts and textbooks on the finite element method e.g. Ref. 22]. In particular 



our notation follows Ref. [22]. Discretizing the differential equation with mentioned Robin 
BCs using the finite element method leads to a system of linear equations. The resulting 
matrix is a banded matrix for which the system of equations can be solved efficiently if high- 
order hierarchical piecewise polynomials are used as a basis and if the degrees of freedom are 
decimated. This hierarchical finite element basis set leads to algebraic systems that are less 
susce pti ble to round-off error accumulation at high order than those produced by a Lagrange 



basis 231]. We use linear hat functions as the linear hierarchical basis. For higher order bases 



we exploit the method of Szabo and Babuska[2J] which relies on Legendre polynomials. 
Below we show the expansion of c{z) in terms of the hat functions and the other higher 
order hierarchical piecewise polynomials on the interval 

p 

C{Z) ^ Q_iiV_i(ei) + CiNii^,) + J2 , (18) 

i=2 
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where = 2{z — Zi)/h + 1; h = Zi — Zi^i and the functions iVj(^) in the interval [—1, 1] are 
given by 

N-Ai) = ^ iV,(0 = ^ (19a) 



^^(0 = l\-,{Odi\ I > 2 (19b) 



These hierarchical bases have useful orthogonality properties that lead to sparse and well- 
conditioned stiffness matrices. Defining an operator C 

C[c] := c"{z) - ^^c{z) (20) 

we can write our differential equation ffTOj) as 

C[c] = r]{z) 

with boundary conditions 

c'{zi) - gc{zi) = ^^^^ 
c'{zu) + gc{zu) = 

The method of weighted residuals is used to construct a variational integral formulation of 
Eq. (12Dl) by multiplying with a test function d{z) and integrating over [zi,Zu]: 

{d, C[c] - T]) = WeH\zi,Zu) (22) 

where is the Sobolev space. We have introduced the inner product 

{d,c) := / d{z)c{z)dz (23) 



Performing the integration by parts in Eq.( !22]) and applying Robin BCs given in Eq.( !2Tl) 
gives 

A{d, c) = {d, 7]) + gd{zi)c{zi) + gd{zu)c{zu) (24) 
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where 

A{d, c) := /" " [-d'{z)c'{z) - -f^d{z)c{z)] dz (25) 

Using the Galerkin approach and exploiting the decimation scheme, we can construct a 
system of hnear equations Be = b where the elements of the vector c are the values of c{z) 
at grid points. The detailed structure of this linear system of equations is given in the 
Appendix A. 



III. NUMERICAL RESULT 



In this section we present the numerical results obtained for the Poisson's solver for 
continuous charge densities with 2DP1DF BC in stand alone mode and for our Ewald-like 
method for point particles interacting by Coulombic potential with 2DP1DF BC. We also 
show numerical evidence for the conservation of energy in molecular dynamics simulation of 
a system composed of sodium chloride atoms. 



A. Numerical results for the Poisson solver 



Our method has an algebraic convergence rate in the non-periodic direction and a faster 
exponential convergence rate in the periodic directions, respectively due to the finite element 
polynomial bases and to the plane wave representation. In Fig. [T] we show the convergence 
rate in non-periodic direction with 7-th order finite elements (p=7 in Eq. flTSl) ). For our 
test, the starting point was the potential rather than the charge density, since the charge 
density can be obtained analytically from the potential by simple differentiation. Our test 
potential had the form 0(r) = sin(asin(^)) sin(6sin(^)) exp(— ^). 



B. Numerical results for point particles 

In this section we give the numerical results of our implementation of the presented 
method for point particles. Since MMM2D is known to be highly accurate we use it as 
reference in this section. First we want to demonstrate that error distribution along the non 



periodic direction is uniform unlike in the 3D periodic methods with correction terms 
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FIG. 1: RMS of relative error for the potential given in Sec. IIII Al with a = 10, b = 10, c = 1. On 
this double logarithmic plot the curve has an asymptotic slope of 14 and machine precision can be 
reached. 
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FIG. 2: Relative error distribution of force norm on each particle along z-axis for 100 random 
systems with 100 particles. 
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171]. To this aim 100 particles were put randomly in a unit cubic cell and the program 



was run 100 times each time with different random positions. Results of the relative error 
of forces exerted on each particle are plotted in Fig. [2l 

In Fig. [3] we show the theoretical scaling 0{N log{N)) can be achieved in practice. The 
crossover with MMM2D for a moderate accuracy of 10~^ in RMS relative error of forces is 
found to be less than 20 particles. Both programs were run in AMD Opteron 2400 MHz. 
The degree of the finite elements is a parameter that can be optimized to obtain the smallest 
possible CPU time for a fixed accuracy. For high accuracies higher degrees are recommended. 
The CPU time for the calculation of the forces dominates in our method over the time needed 
to calculate the energy. 
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FIG. 3: CPU time of one time evaluation of forces on particles and potential energy with our 
metliod(solid curve) and MMM2D method(dashed curve). 

C. Energy Conservation 



Energy conservation is of great importance in molecular dynamics simulations. In order 
to investigate energy conservation in a real simulation, we performed a very long (8 nano 
second) molcular dynamics simulation of a sodium chloride system containing 1000 particles. 
The velocity Verlet algorithm with a time step of 50 atomic units is used to update the 
particle positions and velocities. The short range interactions are obtained from the Born- 



Mayer-Huggins-Fumi-Tosi 25| (BMHFT) rigid-ion potential, with the parameters of Ref. 26 1 
The shortest oscillation period was of the order of 3000 atomic units e.i. 60 molecular 
dynamics steps. After an equilibration for 1 x 10^ steps, 7 x 10® steps were performed during 
which the total energy and potential energy were monitored. The fluctuation of the total 
energy, shown in Fig. HJ has an oscillation amplitude of about 2.5 x 10~^, while the amplitude 
of the potential energy oscillation was 3 orders of magnitude larger. The total energy was 
thus conserved very well. 



IV. CONCLUSION 



In this manuscript we presented a method to solve Poisson's equation for smooth charge 
densities with periodic boundary condition in two directions and finite in the third one. 
It is very efficient for smooth charge densities and it does not require much memory. The 
resulting error distribution is uniform over the entire simulation cell. Our method is based on 
plane wave representation in the periodic directions and finite elements in the non-periodic 
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FIG. 4: The total energy fluctuations calculated with our method. 

direction. Based on this method we can then calculate electrostatic energy and forces of 
particles interacting by Coulombic potential with high accuracy and an A^log(A^) scaling. 
The method satisfies intrinsically and without any approximations the boundary conditions 
approriate for surface problems. It is best suited for a moderate number of particles in 
between 10^ — 10^. The method is expected to be suitable for an efficient parallelization 
since the time dominating parts are only loosely coupled. 
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APPENDIX A: APPENDIX 

We consider a uniform grid on the interval [zi,Zy\ with + 1 nodes {zq, zi, . . . , zn} 
while zq = zi and zn = z^. The interval is thus divided into equally spaced sub inter- 
vals (elements). The functions di^z^ and ci^z) are replaced by the approximate functions -D(-2) 
and Ci^z) which are expanded in the basis of Eqs. ( JT9l) on each subinterval. We use the 
Galerkin approach in which the same bases are used for the expansion of both D{z^ and 
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C{z). Our bases are a combination of the hat function 4''"{z) centered at nodes 



{zj+i - z)/h, z e [zj,Zj+i) 
{z - Zj_i)/h, z e [zj-i,Zj) 
otherwise 



(Al) 
(A2) 



and hierarchical polynomials 2J] (j)'^{z 



otherwise 



(A3) 



localized within the individual elements. Ni are given in canonical coordinates in Eqs.flTQ 
Finally C(z) and -D(-2) within the element \zj_x^Zj\ will be: 



i=2 
P 

D{z) = d,^,<P]_,{z) + d,<P]{z) + J2d,,<P^Az) 



(A4a) 
(A4b) 



i=2 



Note that because (p^ii^) vanishes at all nodes we obtain cj = C{zj). Replacing the 
approximate functions from Eg. ( ]A4al) and Eg. ( ]A4b1) in eguation f l2^ gives 



N 



^[Aj{D, C) - {D,T])j\ = gdoco + gdjyCN 



(A5) 



We split Aj{D,C) as 



A,{D,C)=A^{D,C) + Af{D,C) 



(A6) 
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where 



C{z) within an element is: 



where Cj and 4>j{z) are vectors with p + 1 elements: 



Then 



(A7) 



A^{D,C) := - r D\z)C'{z)dz 

Jzj-l 

Af{D,C) := - r j^D{z)C{z)dz (A8) 
{D,r))j := r D{z)rj{z)dz (A9) 



C{z)=^{z)c, ze[zj.^,zj] (AlO) 



Cj := [cj_i,cj,cj^2,---,cj,pf (All) 



Af(D,C) = cijK.c,- (A13) 
Af{D,C) = d^'MjCj (A14) 



where 

M, ::= - r ^% ^ (A16) 

J Zi-l 

The (p+ 1) X (p + 1) matrix Kj is called the element stiffness matrix and the (p + 1) x (p+ 1) 
matrix Mj is called the element mass matrix. Although the element index j is present in the 
definition of Kj and Mj, in our case of uniform grid spacing these matrices do not depend 

on j. By performing the summation Xlj^i Sj=i ' build up the global mass 

matrix and the global stiffness matrix. We arrange the order of elements of these matrices 
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as: 



Cl 

Co, Cl, . . . , CnY 

Cl,2, • • • , Ci,p, . . . , Cn,2-, ■ 



(A17) 

(A18) 
(A19) 



K 



M 



Kl 
Kc^ 

Ml Mlq 
MIq Mq 



(A20) 
(A21) 



The second term of the summand in Eg. (lASp should be calculated approximately because 
only the values of ri{z) on the nodes are available: 



(A22) 



where 



/,: = 



(f)j{z)r]{z)dz 



(A23) 



Interpolating integration is appropriate to calculate the above integral by fitting a polynomial 
of degree d> 2p to the nodes of element [zj^i, zj] and its neighboring nodes: 



p-i 



(A24) 



Recall that our charge density is localized within the interval [zi,Zu] and it smoothly tends 
to zero at the edges. Therefore it is appropriate to zero pad the ends of the ri{z). The 
coefficients wl are weights from high-order interpolation. Building up the global matrices 
yields: 

{D,r]) = (Fl (A25) 
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where the order of elements of I is the same as in Eq. (lA17p . 



II : -- 



h 

[Jo, /i, . . . , InY 

[h,2-, • • • , h,p-, • • • , In, 2-, • • • , In,i^ 



(A26) 

(A27) 
(A28) 



Finally by adding the right-hand-side of Eq. ( 1A5[) to the global matrices yields: 



Pl 


Mlq 




cl 




II 




^« . 











where Mlq is a sparse (A^ + 1) x N{p — 1) matrix, 

Pq := Kq + Mq 
is a A^(p — 1) X N{p — 1) block-diagonal matrix, 



(A29) 



(A30) 



Pl := -ft^L + Mi geoCQ + geNefj 



(A31) 



is a tridiagonal (A^ -|- 1) x (A^ + 1) matrix, and 



eo := [1,0,. ..,0r 
ejv := [0,...,0,lf 



(A32) 
(A33) 



Multiplying the matrix in Eq. flA29p and eliminating cq in the system of linear equations 
yields: 

[Pl - MlqPq'MIq] cl = Il- MlqPq'Iq (A34) 
Finally we obtain our system of linear equations: 



Bcr = b 



(A35) 
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where the matrix B and the vector b are 

B := Pl- MlqPq'MIq (A36) 
b := h- MlqPq% (A37) 

It turns out that in the general case the matrix B is symmetric tridiagonal of dimension 
(A^ + 1) X [N + 1). The proof for the tridiagonahty of matrix B can be found in the 



context of block cyclic reduction 271]. Note that elements of the vector cl are the values 
of C{z) at the grid points. Therefore by solving a system of linear equations, which has 
a tridiagonal matrix, we can find the values of C{z) at the grid points. Instead of using 
finite element method, we could have used finite differences to solve Eq. (|TOl) . Although 
calculating the right-hand-side b is computationally more expensive in our approach than in 
the finite difference method, the whole process of solving the system of linear equations is 
less expensive because the factorization of the tridiagonal matrix can be done fast. 
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